## Authors: Alexander Herzog and Kenneth Benoit
## Date: May 30, 2015
## Replication file for JOP article "The Most Unkindest Cuts: Speaker Selection and Expressed Government Dissent During Economic Crisis"

rm(list = ls(all = TRUE))

library(arm)
library(lme4)
library(sampleSelection)

### PATHS ########################################
dataAll <- "./generated_data/working_data_all.RData"
dataSub <- "./generated_data/working_data_speakers.RData"
dataOut <- "./estimated_models/all_models_ML.RData"
##################################################

load(dataAll)
load(dataSub)


# Selection models
# ----------------
print("Estimating selection models")    

spoke1.ml <- glmer(data=dataAll, spoke ~ unemployment + safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + crisis + (1 | m), family=binomial("probit"), control=glmerControl(optimizer = "bobyqa"))

spoke2.ml <- glmer(data=dataAll, spoke ~ crisis*unemployment + crisis*safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + (1 | m), family=binomial("probit"), control=glmerControl(optimizer = "bobyqa"))

spoke3.ml <- glmer(data=dataAll, spoke ~ crisis*government*unemployment + crisis*government*safetyScaled + seniorityYearsScaled + party.leader + cabMember + log.party.size + log.debate.days + (1 | m), family=binomial("probit"), control=glmerControl(optimizer = "bobyqa"))


# Outcome models
# --------------
print("Estimating outcome models")    

pos1.ml <- lmer(data=dataSub, textscore ~ unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + (1 | m), control=lmerControl(optimizer = "bobyqa"))

pos2.ml <- lmer(data=dataSub, textscore ~ crisis*unemployment + crisis*safetyScaled + seniorityYearsScaled + government + cabMember + (1 | m), control=lmerControl(optimizer = "bobyqa"))

pos3.ml <- lmer(data=dataSub, textscore ~ crisis*government*unemployment + crisis*government*safetyScaled + seniorityYearsScaled + cabMember + (1 | m), control=lmerControl(optimizer = "bobyqa"))



# Heckman selection model
# -----------------------
print("Estimating Heckman models")    

# ML estimation
joint1.ml <- selection(data=dataAll,
                     selection=spoke ~ unemployment + safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + crisis,
                     outcome=textscore ~ unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis,
                     method="ml"
                     )

joint2.ml <- selection(data=dataAll,
                     selection=spoke ~ crisis*unemployment + crisis*safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days,
                     outcome=textscore ~ crisis*unemployment + crisis*safetyScaled + seniorityYearsScaled + government + cabMember,
                     method="ml"
                     )

joint3.ml <- selection(data=dataAll,
                     selection=spoke ~ crisis*government*unemployment + crisis*government*safetyScaled + seniorityYearsScaled + party.leader + cabMember + log.party.size + log.debate.days,
                     outcome=textscore ~ crisis*government*unemployment + crisis*government*safetyScaled + seniorityYearsScaled + cabMember,
                     method="ml"
                     )

save(spoke1.ml,
     spoke2.ml,
     spoke3.ml,
     pos1.ml,
     pos2.ml,
     pos3.ml,
     joint1.ml,
     joint2.ml,
     joint3.ml, file=dataOut)


